Quantum fluctuations can promote or inhibit glass formation 
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The very nature of glass is somewhat mysterious: while relaxation times in glasses are of 
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sufficient magnitude that large-scale motion on the atomic level is essentially as slow as it is 
in the crystalline state, the structure of glass appears barely different than that of the liquid 
that produced it. 1 ' 2 ' 3 ' 4 Quantum mechanical systems ranging from electron liquids to 
superfluid helium appear to form glasses, but as yet no unifying framework exists 
connecting classical and quantum regimes of vitrification. Here we develop new insights 
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from theory and simulation into the quantum glass transition that surprisingly reveal 
distinct regions where quantum fluctuations can either promote or inhibit glass formation. 

While a wide variety of glassy systems ranging from metallic to colloidal can be 

• • 

accurately described using classical theory, quantum systems ranging from the electronic to 

'X 

magnetic appear to form glassy states. ' Perhaps the most intriguing of these is that the 
coexistence of superfluidity and dynamical arrest, namely the "superglass" state. Recent 
numerical, theoretical and experimental work has given convincing evidence for the reality of 
this unusual arrested state of matter. 7 ' 8 ' 9 However, while such intriguing examples exist there 
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currently exists no unifying framework to treat the interplay between quantum and glassy 
fluctuations in the liquid state. 

To attempt to formulate a theory for a quantum liquid to glass transition we may first 
appeal to the classical case for guidance. Here, a microscopic theory exists in the form of mode- 
coupling theory (MCT), which only requires simple static structural information as input and 
produces a full range of dynamical predictions for time correlation functions associated with 
single particle and collective fluctuations. 10 Although MCT has a propensity to overestimate a 
liquid's tendency to form a glass it has been shown to account for the emergence of the non- 
trivial growing dynamical length scales associated with vitrification. 11 Perhaps more 
importantly, MCT has made numerous non-trivial predictions ranging from logarithmic temporal 
decay of density fluctuations and reentrant dynamics in adhesive colloidal systems to various 
predictions concerning the effect of compositional mixing on glassy behaviour 12 ' 13 . These have 
been confirmed by both simulation and experiment 14 ' 15 ' 16 . 

A fully microscopic quantum version of MCT (QMCT) that requires only the observable 
static structure factor as input may be developed along the same lines as the classical version. 
Indeed, a zero temperature version of such a theory has been developed and successfully 
describes the wave vector dependent dispersion in superfluid helium. 17 In the supplementary 
information, we outline the derivation of a full temperature dependent QMCT. In the limit of 
high temperatures, this theory reduces to the well-established classical MCT, while at zero 
temperature our theory reduces precisely to the aforementioned T=0 quantum theory. The 
structure of these two theories is dramatically different, suggesting the possibility of non-trivial 
emergent physics over the full range of parameters that tune between the classical and quantum 
limits. 
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The fully microscopic QMCT allows for a detailed description of the dynamical phase 
diagram that separates an ergodic fluid region from an arrested glassy one as a function of both 
thermodynamic control variables as well as the parameters (e.g. h ) that control the size of 
quantum fluctuations. To illustrate this, we perform detailed QMCT calculations on a hard- 
sphere system as a function of the system's volume fraction (j) and the dimensionless parameter 
A* (the ratio of the de Broglie thermal wave length to the particle size) that controls the scale of 
quantum behaviour. Despite their simplicity, hard sphere systems are well characterized, 
experimentally realizable, and show all the features of glassy behaviour that are exhibited by 
more complex fluids. It is well-known from experiment and simulation that classical hard 
spheres enter a glassy regime for volume fractions in the range <f> = 50 — 60% independent of 
temperature. 18 ' 19 Figure 1 shows the full structure of the dynamical phase diagram. The QMCT 
calculations are consistent with this in the classical limit (A*— > 0), but upon departure from this 
show a rather remarkable reentrant behaviour. In particular, as the scale of quantum fluctuations 
are tuned from small to high values an initially flat regime is followed by the system becoming 
glassier and then finally favouring the fluid when quantum fluctuations are large. This behaviour 
is surprising given the fact that reentrance is not hinted at in the static structure factor. We also 
show the behaviour produced by a strictly classical MCT calculation performed with the 
quantum structure factor as input where only a featureless border separating liquid from glass is 
demonstrated. This fact clearly shows that the reentrant behaviour predicted by QMCT is a non- 
trivial product of the properties of the theory and not the static structure factor input. 

To obtain a physical understanding of this surprising prediction we turn to the ring 
polymer molecular dynamics (RPMD) approach to quantum dynamics. This method exploits 
the path integral formulation of quantum mechanics in which a quantum particle is mapped onto 

3 



a classical ring polymer consisting of a series of replicas linked by harmonic springs. Static 
properties can be calculated exactly using this mapping while RPMD utilizes the classical 
evolution of the polymers to provide an approximation to quantum dynamics. This approach has 
been previously shown to give accurate dynamical properties for systems ranging from nearly 
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classical to those where tunnelling is dominant. ' 

We performed RPMD simulations for a binary Lennard- Jones system at a density and 
temperature that classically exhibits glassy behaviour (details provided in supplementary 
information). Figure 2(a) shows the change in the diffusion coefficient of the particles as the 
quantum fluctuations of the system, controlled by varying A , are increased. This property shows 
the same reentrance seen in the QMCT results. The structure factor, shown in Figure 2(b), 
reveals only a monotonic broadening over the entire region under study. 

Analysis of the RPMD trajectories allows us to deduce the origin of this effect. In Figure 
2(c), we show the ratio of the average radius of gyration of the polymers representing each 
particle, a static property given exactly by the RPMD simulations, to its free particle value which 
is proportional to A*. The spread of each polymer (or the width of the thermal wave packet) is 
directly related to the quantum mechanical uncertainty about its position. Hence, the uncertainty 
principle dictates that decreasing the width of the packet corresponds to an increase in kinetic 
energy. 

The trend in the spread of the particles as shown in Figure 2(c) is in excellent agreement 
with that seen in the diffusion coefficient (see Figure 2(a)) and provides insight into the reentrant 
behaviour. As quantum fluctuations are introduced into the system the wave packet of each 
particle attempts to delocalize. Initially (A*<0.1) thermally accessible space is available in the 
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system for the particle to expand into, allowing the radius of gyration to increase almost freely. 
The ratio of the spread of the particles to their free particle values is near unity and the diffusion 
is largely unchanged. However as A* is increased further, the width of the packet becomes 
comparable to size of the cage in which it resides. There is now little free space into which the 
packet may expand, leading to a dramatic decline in the ratio. Figure 3(a) shows a configuration 
typical of this regime in which the particle is localized in its cavity. This confinement in its 
position causes a large rise in the kinetic energy of localization exerting, as A* increases, a 
progressively larger quantum pressure on the cavity. 

For diffusion to occur, the particles must rearrange in this highly crowded environment. 
This requires contraction of their wave packets as they pass through the narrow gaps, localizing 
them further and incurring an additional increase in their kinetic energy. This higher energy 
required to push through the gaps acts as a bottleneck to diffusion and leads to the slowing 
reflected in the intermediate A region (see Figure 2(a)). This is also shown in the inset of Figure 
3(a) which depicts the mean square displacement of the particle. A long intermediate beta 
relaxation regime is observed in which the particle is caged before being becoming mobile again 
at long times. When A* is raised further, the thermal wavelength becomes comparable to the 
particle size and the kinetic energy becomes sufficient to flood the barriers between cavities 
leading to a rise in the radius of gyration and the occurrence of tunnelling between the cavities, 
thereby facilitating diffusion. This can be seen in the representative snapshot shown in Figure 
3(b) in which the particle is stretched across two cavities. Accordingly, the ratio of the radius of 
gyration to its free value recovers with a corresponding increase in diffusion and diminishing of 
the caging regime as shown in the inset of Figure 3(b). 
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The theory we have developed for the quantum glass transition predicts interesting 
generic dynamical anomalies such as a reentrant border between the disordered arrested and fluid 
regimes. Semi-classical quantum dynamics simulations display similar features and physically 
illuminate the origin of the predicted relaxation motifs. The physical interplay between crowding 
and quantum derealization reported here, a generic feature of quantum glassy systems, might 
also be responsible for other physical phenomena. For example, it has been experimentally 
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observed that lighter isotopes of hydrogen diffuse more slowly than heavier ones in water and 
palladium 25 which has been recently elucidated by theory. 26 The regime where reentrance is 
observed is therefore realizable in chemical systems near ambient temperatures. 

It is likely the reentrant transition observed here may also have implications beyond 
glassy systems. Intuition suggests that increasing quantum fluctuations monotonically enhances 
the exploration of the energy landscape. This forms the basis of the quantum annealing approach 
to optimization ' . However, our work indicates that in certain regimes increasing quantum 
fluctuations can lead to dynamical arrest and hinder optimization. Indeed, reentrance has recently 
been observed in the dynamical phase behaviour of simple models of quantum optimization 
currently under investigation in the field of quantum information science. 29 Thus deep 
connections exist that unite these seemingly distinct physical systems and processes. 
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Figure 1: Dynamic phase diagram calculated from the QMCT for a hard-sphere fluid. <f>is 



the volume fraction, A* = ^K l jma'- is the thermal wavelength in units of inter-particle 
separation a , and (5 = l/k B T is the inverse temperature. Inset: Dynamic phase diagram using the 
quantum mechanical input with a classical MCT. 
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Figure 2: Diffusion as function of quantumness from RPMD simulations, (a) The diffusion 
constant as a function of the quantumness, A* obtained from the RPMD simulations for a 
quantum Kob-Anderson Lennard Jones binary mixture for T* = 2.0 (red curve) and T =0.7 
(black curve) . (b) Classical and quantum static structure factor of the "A" type particles, (c) 
Root-mean-square of the radius of gyration as a function of A =0 for the two systems shown in 
panel (a). The radius of gyration is defined as the average distance of the replicas from the 
polymer centre. 
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Figure 3: Snapshots from the RPMD simulations taken from the caged (A =0.7, panel (a)) 
and tunnelling (A*=1.3, panel (b)) regimes at T* = 0.7. For clarity all but one ring polymer in 
each snapshot is represented by its centre of mass. The red spheres represent the replicas of the 



polymer. Insets depict the mean square displacement, ( R(t) - R(0) 



calculated from RPMD 



(solid curves) and classical MD (dashed curves) in Lennard Jones reduced units. The inset of 
panel (a) shows a long beta relaxation regime compared to the classical simulation and the 
tunnelling case shown in (b). 
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The quantum mode coupling theory 

We first outline the derivation and explicit expressions of the quantum mode coupling theory 
(QMCT). As discussed in the caption of Fig.l, the magnitude of quantum fluctuations may be 
measured by a dimensionless parameter that sets the ratio of the thermal wave length, A to the 
particle size, a namely: 



where H is Planck's constant divided by 2tt, k B is Boltzmann's constant, T is the temperature, m is 
the particle mass, and a is the particle diameter. 

Note that even for hard-spheres a temperature appears that defines the scale of kinetic fluctu- 




A 



a 



(1) 
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ations. This may be thought of as an intrinsic noise temperature following the quantum fluctuation- 
dissipation theorem (QFDT). Defining a projection operator based on the Kubo transformation, 

(A\B) = —J d\(A(-i\)B(0)) = — J dX(A(0)B(iX)), (2) 

where time evolution is defined via the standard Heisenberg picture: 

A(t) = e iHt l h Ae- iHt/r \ (3) 

an exact quantum mechanical equation of motion is found for the Kubo-transform of the density- 
density correlation function, 

4> q (t) + &Mt) + I dt'M q (t')^ q (t - t') = 0, (4) 
Jo 

where D? q = J^^ g 2 , <j> q (0) is the Kubo-transformed static structure factor, and M q (t) is the 
memory function. 

At high temperature the frequency coefficient in the second term reduces to 



l k B Tq 2 
/-x""' \\ mS, 



lim tt q = a — , (5) 



and at low temperature becomes, 



hq 2 

^to q = -—-=u q (6) 



r^o q 2mS q 



which may be recognized as the well-known Bijl-Feynman dispersion at zero-temperature.GE 



Following the mode-coupling approach generalized to the quantum mechanical context, the 
following expression for the memory function of Eq. Q may be derived (in what follows we use 
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the notation C(cu) = dte llJJt C{t) for quantities in frequency space): 

*.<«> « £> u < 7 > 



k 

x(u> - u/)T(u/, u - u)')(j) q ^ k {u')(j) k {u - u'), 



where 



T(ui,u 2 ) = n(-wi)n(-w 2 ) - n(a;i)n(o;2), (8) 
and the vertex is given by 



v q (k,q- k) 



An(n q _ k )An(n k )C, 



q,k,q—k 



Sq-kSkK(Q q _ k , fl k ) 

with 



(9) 



r , n q S q S k S q - k - [q ■ kSq. k + q-{q- k)S k ] 

q ' k ' q ~ k " Q q An(Q k + n q ^ k )-(Q k + Q^ k )An(n q ) 

Here K(Q q _ k , Q k ) = + An( W ) = n(w) - n(-w) and n( w ) = ^ is the 

Bose distribution function at temperature T. 

The above expressions close the equation of motion [4], and require only the static structure 
factor to produce a full approximation to the time dependence of the quantum density-density time 
autocorrelation function. To derive the full expressions quoted above we have resorted to a finite 
temperature generalization of the "resonance approximation" of Gotze and Liicke. 3 4 However, the 
results presented in Fig. 1 (and thus the predicted reentrance effect) are robust and do not depend on 
this approximation. This can be demonstrated by substituting a variety of different approximations 
to remove the dependence of various static terms in the verticies on the integrations over imaginary 
time induced by Kubo transformation. 
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Using the input of accurate quantum structure factors (as described in the next subsection) 
one can make predictions as to the role of quantum fluctuations on the glass transition. A version 
of QMCT has previously been developed to treat the quantum liquid regimeP This theory is not 
capable of treating the regime where dynamics become glassy. A future article will detail the 
relationship between the theory used here and the previous version of QMCT. 

It may be shown analytically that the above equations reduce to the venerable classical mode- 
coupling equations in the high temperature limit and to the Gotze Lucke theory at T = 0. The latter 
theory produces a representation of the dispersion of superfluid helium that is at least as accurate 
as the Feynman-Cohen (FC) theory ® at low values of q and exhibits Pitaevskii-bending of the 
spectrum at high q, unlike the FC theory. In particular at high T, 

\im M q (t) = ^ n - [ £k(q-kc k + q-(q-k)c q „ k f<P q - k (t)<p k (t) } (11) 

where n is the number density and c q = ^{l — J-l is the direct correlation function. In addition, 
4> q (t) reduces to the classical intermediate scattering function, F(q, t) as f3 — > 0. This is recognized 
as the classical MCT memory function.^ 

At T = the equation for the memory function reduces to: 

Mq{u) = £j / Wf V " {K q ~ k) Io " - d2) 

with 

Ti 

v q (k, q - k) = 7^{u k + u q - k + u) q ) {q-kc k + q-{q- k)c q ^ k ) (13) 

which are the T = equations for quantum density fluctuations in superfluid helium first derived 
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by Gotze and Liicke. 3 4 Note that in the T = case, the entire structure of the memory function 
differs greatly from that of its high temperature counterpart and the convolution structure is lost. 



Eqs.(|T2[) and ( [13] ) do not imply a memory function that is a product of correlators at identical 
times. This is a consequence of the QFDT that must be satisfied. At T = the function T(u q , u>k) 
becomes proportional to the difference of a product of step-functions in frequency, dramatically 
altering structure of the theory. This distinction between the low and high temperature limits has 
important consequences. In addition to the robust prediction of reentrance, we also find that glassy 
behavior cannot be supported in the strict T = case. Some or all of these features seem to emerge 
both in certain quantum spin glasses and in recent work on quantum versions of lattice models of 
glassy liquids 8 . A future paper will be devoted to both a more explicit derivation of the theory 
outlined here as well as the physical implications of our work and the connection to other models 
of quantum glass behavior. 



Quantum integral equations for static structure 

The quantum integral equation approach used in this work to generate the input required by the 
QMCT is based on the early work Chandler and RichardsonJ^H^I For completeness, we provide an 
outline of the approach. We begin with the Ornstein-Zernike relation applicable to the quantum 
liquid. The quantum system composed of N particles can be mapped on a classical system consist- 
ing of N ring polymers, each polymer being composed of P beads. Then, we can write the matrix 
RISM (reference interaction site model equation for the classical isomorphic system by: 

h(\r — r'|) = u * c * co(\r — r'|) + ntu * c * h(\r — r'|), (14) 
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where * denotes a convolution integral and n is the number density. In the above equation, h(r), 
u)(r), and c(r) are the total correlation function, the self correlation function, and direct correlation 
function, respectively, defined by: 

1 fW 
h(r) = — j dXh(r, A) 

i r m 

uj(r) = — / dXu(r,X) (15) 



1 



1,3 



c(r) = ^ / dAc(r, A) 



o 



and /i(r, A), u;(r, A), and c(r, A) are the imaginary time total, self, and direct correlation functions, 



respectively. In the classical limit Eq. ( 14) reduces to the classical Ornstein-Zernike equation with 
u(r) = 1. In what follows, we will use the notation cu 5 (A) for the Fourier transform of co(r, A), and 
similarly for c q (X) and h q (X): 



1 



tiB 



h <i = jrp I d\h q (X) 











hd 

dXu q (X) (16) 



1 

hp 

We now use the mean-pair interaction approximations along with the quadratic reference 
action 9 and rewrite: 

Co q (X)=exp{-q 2 R 2 (X)} } (17) 

where 

^l-co^A) 

v ' ^ (3mQ 2 + a- 
j j j 

m is the particle mass, Vtj = 2nj /h(3 is the Matsubara frequency and <x, is given by: 

I poo i-hfi 

ocj = dq dXq\{l - cos{QjX)Cj(q,X). (19) 

07T Hp Jq Jq 
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In the above the solvent induced self-interaction is given by: 

Vq = -c 2 q {noj q + n 2 h q ). (20) 

We now need to close the quantum Ornstein-Zernike equations, which in g-space can be 
written as: 

h g = UqCgUq + nu q c q h q . (21) 
We use the Percus-Yevick closure of the form (in r-space): 

c(r) = (h(r) + c(r) + l)(exp(-fiv(r)) - 1), (22) 

where v(r) is the pair interaction between two particles. 

RPMD Simulations 

We performed RPMD simulations of the Kob-Anderserf^ binary Lennard- Jones (LJ) glass forming 
system. The Lennard Jones potential between particles i and j is given by, 



(23) 



The parameters and their conversion to atomic units as used in this work is given in Table [T] The 
systems consisted of 1000 particles, 800 of type A and 200 of type B in a cubic box of length 9.4 
oaa- The equations of motion were integrated using a timestep of 0.005 LJ units (0.35 fs) using the 
scheme of reference^. The simulations were carried out at constant volume for consistency with 
the QMCT results. The RPMD simulations were performed for distinguishable particles which is 
a valid approximation in the regime where the reentrance is observed. 
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The number of beads, P, used was given by the formula, 



P = ^ (24) 



which was found to give good convergence for all the regimes studied. 

Initial configurations were generated by annealing from a temperature T*=5.0 to the target 
temperature over a period of 1,000,000 timesteps. From these initial configurations we ran a further 
200,000 steps of equilibration using the a targeted Langevin equation normal mode fhermo starting 
schemePl This was followed by microcanonical dynamics for 2,000,000 steps during which the 
results were collected. The quantum effect, A*, was varied by changing the parameter h. Five 
simulations were run for each temperature and value of h and the results averaged. 

The root mean square radius of gyration is defined as, 

P V V2 



' P 

k=l 



pEltf'-'Pl') • <25) 



where 

p 



r^IM" (26) 

fe=l 



is the center of mass of the ring polymer representing particle i. In the free limit the radius of 
gyration is 



r 1 / h2 

r?-2^r x > 2 <27) 

which is related to the De Broglie thermal wavelength as defined in Eq. [T] via multiplication by 
two. 
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Parameter LJ units Atomic Units 

e AA 1 3.8xl0" 4 

e BB 0.5 1.9xl0~ 4 

e AB 1.5 5.7xl0~ 4 

a AA 1 6.43 

a BB 0.88 5.65 

u AB 0.8 5.14 

Mass A 1 3646 

Mass B 1 3646 



Table 1: Parameters used in our RPMD simulations on the Andersen-Kob Lennard- Jones glass 
forming system. 
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